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ABSTRACT 

Coherent oscillations have been observed during thermonuclear (Type I) X-ray bursts from 14 accreting neu- 
tron stars in low mass X-ray binaries, providing important information about their spin frequencies. However, 
the origin of the brightness asymmetry on the neutron star surface producing these oscillations is still not un- 
derstood. We study the stability of a zonal shearing flow on the neutron star surface using a shallow water 
model. We show that differential rotation of > 2% between pole and equator, with the equator spinning faster 
than the poles, is unstable to hydrodynamic shear instabilities. The unstable eigenmodes have low azimuthal 
wavenumber m, wave speeds 1 or 2% below the equatorial spin rate, and e-folding times < 1 s. Instability 
is related to low frequency buoyantly driven r-modes that have a mode frequency within the range of rotation 
frequencies in the differentially rotating shell. In addition, a modified Rayleigh's criterion based on potential 
vorticity rather than vorticity must be satisfied. The eigenfunctions of the unstable modes have largest ampli- 
tudes nearer to the rotational pole than the equator. We discuss the implications for burst oscillations. Growth 
of shear instabilities may explain the brightness asymmetry in the tail of X-ray bursts. However, some fine 
tuning of the level of differential rotation and a spin frequency near 300 Hz are required in order for the fastest 
growing mode to have m = 1 . If shear instabilities are to operate during a burst, temperature contrasts of > 30% 
across the star must be created during ignition and spreading of the flash. The frequency drift may be related 
to the non-linear development of the instability as well as the cooling of the burning layer. The properties and 
observability of the unstable modes vary strongly with neutron star rotation frequency. 
Subject headings: accretion, accretion disks — X-rays:bursts — stars:neutron 



1. INTRODUCTION 

Type I X-ray bursts are thermonuclear flashes on the sur- 
face of accreting neutron stars in low mass X-ray bina- 
ries (LMXBs) (see Lewin, van Paradijs, & Taam 1993, 
1995; Strohmayer & Bildsten 2003 for reviews). Nearly- 
coherent oscillations with frequencies in the range 45-620 Hz 
have been seen during Type I X-ray bursts from 14 sources 
(Strohmayer et al. 1996; Strohmayer & Bildsten 2003 and 
references therein; Strohmayer et al. 2003; Kaaret et al. 2003; 
Villarreal & Strohmayer 2004). It is thought that rotational 
modulation of a brightness asymmetry on the neutron star 
surface produces an oscillation in X-ray intensity at the spin 
frequency of the neutron star (e.g. Strohmayer & Markwardt 
1999). This has provided spin measurements for these LMXB 
neutron stars, confirming their rapid spin, as expected if they 
are the progenitors of the millisecond radio pulsars (Bhat- 
tacharya 1995). 

The observational properties of burst oscillations have 
been intensively studied. The oscillations are sinusoidal 
(Strohmayer & Markwardt 1999; Muno et al. 2000); they drift 
upwards in frequency by a few Hz during the burst, reaching 
an asymptotic frequency that is stable to a few tenths of a per- 
cent on timescales of years (Strohmayer et al. 1998b; Muno 
et al. 2002a); they are observed when the source is accreting 
at higher accretion rates rather than lower (Muno et al. 2000; 
van Straaten et al. 2001; Franco 2001; Muno, Galloway, & 
Chakrabarty 2004); they are seen in the burst tail as well 
as the rise, with amplitudes typically 5% in the tail (Muno, 
Ozel, & Chakrabarty 2002) and much larger during the burst 
rise (Strohmayer, Zhang, & Swank 1997b; Strohmayer et 
al. 1998a); their amplitude grows with photon energy (Muno, 
Ozel, & Chakrabarty 2003). That the oscillations occur at 
the stellar spin frequency has been confirmed by detection of 
burst oscillations from two of the accreting ms X-ray pulsars 



(Chakrabarty et al. 2003; Strohmayer et al. 2003). Burst oscil- 
lations were also seen during the hours-long superburst from 
4U 1636-54 (Strohmayer & Markwardt 2002). 

The physics of the burst oscillations has, however, re- 
mained unsolved. Initial work concentrated on the upwards 
frequency drift, which Strohmayer et al. (1997a) suggested 
was related to angular momentum conservation in the burn- 
ing layer, which expands outwards when heated and slowly 
contracts during the cooling tail of the burst. However, de- 
tailed models of the vertical expansion (Cumming & Bild- 
sten 2000; Cumming et al. 2002) showed that the expected 
spin down from this effect was a factor of 2-3 smaller than 
observed. Spitkovsky et al. (2002) studied the ignition and 
propagation of a hotspot during the burst rise. They argued 
that rapid rotation is crucial because it allows the horizon- 
tal pressure gradients associated with a localized hot spot to 
be supported by hurricane-like winds. They showed that this 
leads to a new mode of burning front propagation in which 
the ageostrophic fluid flow associated with expansion of the 
hotspot brings in fresh fuel, allowing the burning to propa- 
gate, and spread around the star on a timescale of a second or 
less. 

The picture of a spreading hotspot during the burst rise 
is appealing, and consistent with the large amplitudes ob- 
served (Strohmayer et al. 1997b). It does not however explain 
why burst oscillations are observed in the cooling tail, when 
the whole surface has ignited. There have been two mech- 
anisms suggested that might cause asymmetries in the burst 
tail. Spitkovsky et al. (2002) pointed out that the spreading 
burning front leaves behind temperature variations across the 
surface of the star, with associated zonal flows. They specu- 
lated that these flows were unstable to vortex formation, and 
pictured a vortex trapped in the zonal flow somewhat analo- 
gous to Jupiter's Great Red Spot. 
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The second suggestion is that oscillation modes are excited 
in the burning layer, and that the frequency drift is due to a 
changing mode frequency as the burning layer cools (Heyl 
2004; Lee 2004; Piro & Bildsten 2005; see also McDermott & 
Taam 1987; Strohmayer & Lee 1996; Piro & Bildsten 2004). 
In particular, Heyl (2004) identified low azimuthal wavenum- 
ber, buoyancy-driven r-modes as promising candidates, since 
they drift backwards in the rotating frame, and have a rel- 
atively large latitude coverage. However, the frequency of 
the lowest radial order mode in the burning layer is too large 
to explain the observed drifts, and so there remained a ques- 
tion as to whether a mode with the correct frequency exists, 
and how it is excited. Piro & Bildsten (2005) answer the first 
question by proposing that the surface wave associated with 
the burning layer transitions to a crustal interface mode dur- 
ing the burst tail, by an avoided crossing as the layer cools. 
Because the crustal interface mode has a frequency that is in- 
dependent of the temperature in the burning layer, this termi- 
nates the frequency drift at a frequency that depends only on 
properties of the crust/ocean interface. 

In this paper, we connect the ideas of unstable zonal flows 
and oscillation modes by studying the linear stability of shear- 
ing zonal flows during the tail of a Type I X-ray burst. We 
show that for differential rotation between pole and equator 
of > 2%, there are linearly unstable modes with properties 
that may be relevant for understanding burst oscillations. The 
properties of the unstable modes are closely related to the r- 
modes of a rigidly rotating layer. The unstable modes have 
low azimuthal wavenumbers m, growth rates < 1 s, and wave 
speeds of order 1% smaller than the equatorial rotation rate. 
We start in §2 by discussing the zonal flows expected during 
the tail of a burst, and then in §3 describe calculations of the 
linear stability of a quasi-geostrophic shallow water model of 
the burning layer. We then discuss the implications for our 
understanding of burst oscillations in §4. In the Appendix, so- 
lutions to the full set of shallow water equations are presented, 
and compared with the quasi-geostrophic solutions calculated 
in §3. 

2. ZONAL FLOWS DURING THE TAIL OF A BURST 

Spitkovsky et al. (2002) emphasized that a zonal flow natu- 
rally arises as the burning front spreads around the star. If the 
local cooling time is the same at each point, the part of the sur- 
face that ignited first will be cooler by a factor of ~ t sprea d /? coo i 
than the part that ignited last, where f spre ad is the time for the 
burning front to spread across the surface, and f coo i is the lo- 
cal cooling time. Because the spin period of the star is much 
shorter than either of these timescales, Coriolis forces act to 
produce a zonal flow which balances the lateral pressure gra- 
dient associated with the differential cooling. 

We start by estimating the amount of differential rotation 
expected. We model the burning layer using the shallow water 
equations (e.g. Houghton 1986; Pedlosky 1987), in which the 
detailed vertical structure of the layer is neglected. We work 
in the rotating frame. The momentum equation is 
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where u = egug(9, (\j) + e$u^(9 ', (j>) is the horizontal fluid veloc- 
ity, H is the pressure scale height of the layer, g = GM/R 2 is 



the surface gravity and fl is the neutron star spin frequency. 
Throughout this paper, 9 is chosen such that 6 = corresponds 
to the pole, and 9 = 90° the equator. The advective deriva- 
tive d/dt = d/dt + v'V involves horizontal velocities only, 
and also only Coriolis forces arising from horizontal motions 
are included in equation {0- In addition, we assume that the 
burning layer is spherical, and neglect distortions due to the 
centrifugal force. 

We assume a steady state flow with azimuthal velocity u = 
u,pe<b, where we write UA9) = Ail(8)rsin9. The force balance 



is 



2flcos9u c i 



or in terms of /i = cos 9, 
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We will consider in this paper the simplest shear profile sym 
metric about the equator, 

in which case integrating equation @ gives 
H(ji) = H (l + ¥-n 4 
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where Ho is the scale height at the equator. The dimensionless 
parameter F is 
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where R = 10 R(, km and g = g\4 10 14 cm s" 2 . Another way 
to write this is F = (R/Rd) 2 , where in terms of the Brunt- 
Vaisala frequency N, Rd ~ H(N/fl) ~ 1 km is the Rossby 
deformation radius, the horizontal scale on which rotational 
effects are important (Pedlosky 1987). 

In equation Q, we write the pressure scale height H in 
units of 10 3 cm. At the high temperatures achieved during 
the burst, the ideal gas equation of state is appropriate, and ra- 
diation pressure may also play a role at the peak temperature 
(Cumming & Bildsten 2000). Therefore, we estimate 
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where \L m is the mean molecular weight (we use a subscript 
m to avoid confusion with /i = cos 9), and f3 is the ratio of gas 
to total pressure. For typical fluxes and base temperatures ap- 
propriate for the peak of a burst, Cumming & Bildsten found 
(3 w 0.5, with (3 rapidly approaching unity as the layer cools. 
Then, 

F=m ta£(_^_\ R 2 (9) 
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where is the temperature in units of 10 9 K. 

Rewriting equation (j5Jl, the differential rotation between the 
pole and equator is 



(10) 



where 0| = GM/R 3 = g/R (this is a similar result to 
Spitkovsky et al. 2002, eq. [74]). This gives 
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If s > 0, the equator is rotating faster than the poles, and the 
fluid piles up at the poles, AH > 0. 

If we follow Spitkovsky et al. (2002) and take AH/H ps 
fspread/W, we find s w 1% (f spre ad/l s)(10 s/f cool ). The direc- 
tion of the differential rotation will depend on the details of 
ignition and spreading of the burning. Spitkovsky et al. (2002) 
argue that ignition starts at the equator, in which case at any 
given time as the burst cools the poles will be higher pressure 
relative to the equator. This implies s > 0, or more rapid ro- 
tation at the equator than the poles. In fact, this case is more 
favorable for shear instabilities. We will see in §3 that insta- 
bility sets in for much smaller levels of differential rotation if 
the equator is rotating more rapidly than the poles. 

3. LATERAL SHEAR INSTABILITIES IN THE SHALLOW 
WATER MODEL 

3.1. Perturbation equations with the geostrophic 
approximation 

We now perturb the shallow water equations around the 
background zonal flow described in ^2] To simplify the cal- 
culation, and to highlight the physics of the instability, we 
will assume that the perturbed fluid velocity is in geostrophic 
balance. This approximation must in fact break down at the 
equator, where the vertical component of the rotation vector 
vanishes. However, away from the equator, we expect the ap- 
proximation to be a good one. In the Appendix, we calculate 
the r-modes of the full set of shallow water equations, and 
show that they compare well with the geostrophic solutions 
except for a small region near the equator. 

The assumption of geostrophic balance for the perturba- 
tions allows us to work with the vorticity equation only. Tak- 
ing the curl of equation Q, and using equation (0, gives the 
vorticity equation for shallow water, 



dt { H 



= 0, 
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where u r = ( Vx u) r is the radial component of fluid vorticity 
in the rotating frame, 2£l/i is the radial component of vortic- 
ity coming from the overall rotation, and H is the thickness 
of the layer. For H constant, this is simply a statement of 
conservation of vorticity; with a variable H, the potential vor- 
ticity (u) r +2flp)/H is conserved. The extra physics included 
in potential vorticity conservation (rather than ordinary vor- 
ticity conservation) is vortex stretching: as a column of fluid 
is stretched vertically, its vorticity changes to satisfy overall 
angular momentum conservation (e.g., Houghton 1986; Ped- 
losky 1987). 
The vorticity of the background state is 

n r = -(d/d^)[(n+An)(\-^ 2 )\ , (13) 

which is Q, r = 20yu for rigid rotation, and Q, r = 2f2/i(l +s) - 
4fX?/i 3 for the rotation law of equation Perturbing the 
vorticity equation then gives to linear order in the perturba- 
tions 
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where H and Vl r are the height and vorticity of the back- 
ground, and h and ui r are the perturbations to the height 
and vorticity, and the perturbation to the fluid velocity is 

Su^e^ + Sugeg. 



Assuming that the perturbations are in geostrophic balance 
provides a relation between the fluid velocities and h. We first 
write the perturbed velocity in terms of a stream function tjj, 
such that 
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which automatically enforces V-u = 0. Then uj r = -Cip/r 2 , 
where 
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is the Legendre operator, and m is the azimuthal wavenumber, 
i.e. we have assumed a <fi dependence oc e" n ^ . We then use the 
cp component of the momentum equations under the assump- 
tion of geostrophic balance, 2flcos9ug = (-g/ r sin 9)(dh/d<fi), 
to derive a relation between ip and h, 
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It is important to note that this relation is not exact, because 
it only approximately satifies the 9 component of the momen- 
tum equations. The approximation is that geostrophic balance 
holds at the equator, which it does not. We also note here that 
this approximation removes gravity waves from the system. 

We now look for solutions ip oc exp(-imuit + imcj)), where 
Q = u> + ia is the complex mode frequency. Also, we write the 
perturbation equation in dimensionless variables, with unit of 
time Or 1 = 0.53 ms (/i/300 Hz)" 1 , and velocity flR. The final 
perturbation equation is 



(uj-AU)£ip-il>— — — = 0, 



(18) 



with C given by equation fl!6l >. With appropriate boundary 
conditions, this specifies an eigenvalue problem for the com- 
plex frequency u>. 

We solve the eigenvalue problem with a shooting method. 
We reduce equation dl 81 to two first order differential equa- 
tions for ip and $ = (1 - /j, 2 )(dip/dij). Then, following Dikpati 
& Gilman (1999) and Bildsten, Ushomirsky, & Cutler (1996), 
we remove the singularity at p, = 1 by defining / and g such 
that V = (1 -Ai 2 ) |m|/2 / and $ = (1 - n 2 ) l "' l/2 g, which also en- 
forces the boundary condition ^ = at the pole. The resulting 
equations are 



(iV)^ = s+HH/ 

dfx 



and 
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We start integrating a small distance from the pole, by setting 
/=lat^t=l-e and g = — p, \m\ f. For our choice of rotation 
law, the equations are symmetric under the change /i — > — fi. 
Given this, we integrate to the equator (p, = 0) and look for 
solutions which are either symmetric or antisymmetric about 
the equator (either g = or / = at fi = 0). Since the eigen- 
value ui + ia is complex, we integrate both real and imaginary 
parts of / and g, and perform a two-dimensional search over 
us and a using Newton's method. 
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FIG. 1. — Wave speeds and growth rates as a function of the differential 
rotation s for 2D modes (F -C 1). The background vorticity profile satisfies 
Rayleigh's condition for instability for s > 0.2, or to the right of the dotted 
line. The dashed lines demark the "corotation region" in which a critical 
level exists in the flow. The two solid lines are for modes which are the 
/ = 2, m = 1 and / = 4, m = 3 r-modes at zero differential rotation. The 1 = 4 
mode approaches the boundary of the corotation region and stops. The 1 = 2 
mode enters the corotation region and becomes unstable. The long-dashed 
curve shows the growth rate (multiplied by a factor of 10). 



3.2. Shear instabilities for a two-dimensional flow 

We present the results of the stability calculation in 33.31 
First, it is useful to summarize the results of stability of lat- 
itudinal shear in a two-dimensional spherical shell. In the 
shallow water model, the limit of a 2D shell is equivalent to 
holding H constant, or setting F = 0. This case has been stud- 
ied by Watson (1981), Dziembowski & Kosovichev (1987), 
Charbonneau, Dikpati, & Gilman (1999), Watts et al. (2003), 
Watts, Andersson, & Williams (2004), and in the weakly non- 
linear regime by Garaud (2001). Much of this previous work 
was motivated by trying to understand the rotation profile of 
the Sun, in particular the role of shear instabilities in the so- 
lar tachocline. These results can be understood in a simple 
way which is helpful when discussing the shallow water in- 
stability. Dynamical instability is closely connected with the 
existence of a critical level at which the real part of the mode 
frequency is equal to the local rotation rate, as found for shear 
instabilities in differentially rotating disks (e.g. Papaloizou & 
Pringle 1985) (also, the existence of a critical level leads to a 
complex mode spectrum, including a continuous spectrum of 
neutral modes, see e.g. Watts et al. 2003). 

3.2. 1 . Rayleigh 's criterion 

Watson (1981) derived necessary conditions for instabil- 
ity starting with equation Jl 8I > with F = 0. Multiplying by 
^*/(<D-Af2), integrating over fi, and taking the imaginary 
part leads to the requirement that dn r /d/j, must change sign 
somewhere in the domain for instability. This is the equiv- 
alent of Rayleigh's criterion for instability of parallel shear 
flows (e.g. Drazin & Reid 1981), that there be an inflexion 
point in the velocity profile 1 . For the rotation law of equa- 
tion l|5}, dQr/dfi = 2(1 +s-6s/j, 2 ), which gives the necessary 

1 Note that there is an additional necessary condition for instability, 
Fj0rtoft's criterion, that is always satisfied for our choice of rotation law. 



FIG. 2. — The critical value of s needed to satisfy Rayleigh's criterion for 
the shallow water equations as a function of F. Solid line is for the rotation 
law fi = 1-i/i 2 , dotted line is for fi = l-(s/2)(> 2 +/i 4 ). For F < 1, s crit = 1/5 
and 1/7 respectively for these rotation laws. For F>1, j cr j t fa 8/3F and 
48/19F respectively. 



condition for instability s>l/5ors<— 1. Note that insta- 
bility is possible for much smaller differential rotation when 
the equator is rotating faster than the pole, as opposed to the 
pole rotating faster than the equator. For this reason, we will 
concentrate on the s > case in this paper. 

It has been pointed out in studies of the Sun's rotation pro- 
file that introducing a /i 4 term in the rotation law leads to in- 
stability at a lower level of differential rotation than with a /i 2 
term only (Charbonneau, Dikpati, & Gilman 1999). This can 
be understood from Rayleigh's criterion. For example, the ro- 
tation law £1= l-(s/2)(/i 2 + /i 4 ) considered by Garaud (200 1 ) 
requires s > 1/7 = 0.14 for Rayleigh's criterion to be satisfied, 
rather than s > 1 /5 as previously, allowing instability to set in 
at a lower differential rotation. Interestingly, this is not the 
case for s < 0, for which the requirement is still s < — 1 even 
when a [f term is included. 

3.2.2. The r-modes, critical levels and instability 

A rigidly-rotating 2D shell supports a set of r-modes, with 
frequency in the rotating frame 

2mfl 

(21) 



/(/+!)' 



and eigenfunctions P'"(fi), with m < I (to see this, set F = 0, 
Ail = 0, and dil, /dfi = 20, in eq. 1181 ). The inertial frame 
frequency is m(il+ uS). For m=l, the modes with / odd (even) 
are symmetric (antisymmetric) about the equator, and oppo- 
site for m = 2. As described by Watts et al. (2003), the in- 
stability can be understood in terms of the evolution of these 
modes as the level of differential rotation is increased. In Fig- 
ure[flwe show two examples of the evolution of the mode fre- 
quency with increasing s, focusing on the s > case. The crit- 
ical value for Rayleigh's criterion s = 1/5 is indicated by the 
dotted line. As the level of differential rotation is increased, 
the onset of instability is associated with the development of 
a "critical level" in the flow, at which the mode frequency is 
equal to the local rotation frequency, or li = il(/i). This oc- 
curs for points within the wedge defined by the lines ui = 1 
and uj = l-i, shown as short-dashed lines in Figure 
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FIG. 3. — The r-mode frequencies for m = 1 (solid) and m = 2 (dashed) 
modes with I - 2, 3,4, 5 as a function of F for rigid rotation (s = 0). For 
large F, the scaling is ui oc F~ 1 / 2 , as p redicted by the asymptotic analysis of 
Longuet-Higgins (1968) (see eq. |Z5I ), The dotted curve shows the critical 
value of .s needed for Rayleigh's criterion to be satisified at a given value of 
F. Only modes which lie above the dotted line are expected to be unstable. 
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FIG. 4. — Antisymmetric m = 1 r-mode eigenfunctions, for s = 0, F = 100. 
We show the height perturbation h in one hemisphere, normalized to be unity 
at the maximum value. These modes have u = 1,3,5 and 7. The frequencies 
are -0.0682,-0.0312,-0.0211, and -0.0160. The predicted turning points 
for these modes are at fi = (2u+ l) 1 / 2 //" -1 / 4 (see Appendix) which gives 57° , 
and 33° for v — 1 and 3, in good agreement with the numerical solutions. The 
remaining two modes are able to propagate all the way to the pole. 



The unstable modes are those which, as differential rotation 
increases, satisfy Rayleigh's criterion at the point at which 
they first develop a critical level (Watts et al. 2003). The solid 
curves in Figure [2 show the evolution of the I = 4, m = 3 and 
l = 2,m=\ r-mode frequencies as differential rotation is intro- 
duced. The / = 4 mode approaches the corotation region when 
the background vorticity profile has not yet met Rayleigh's 
criterion (to the left of the dotted line), and does not cross 
the corotation boundary. The I = 2 mode enters the corota- 
tion region with Rayleigh's criterion satisfied (to the right of 
the dotted line), and becomes an unstable mode. The long- 
dashed curve in Figure [2 shows the growth rate of the mode, 
which increases with increasing s. 

We find instability for s > 0.29, in agreement with Watson 
(1981). The reason for this is that the unstable modes have 
uj w 0.7 because they are related to the 1 = 2 r-mode, and so 
s > 0.3 is needed in order for a corotation point to exist in 
the flow. This differential rotation profile has two modes that 
are unstable, the I = 2,m = 1 antisymmetric mode, and the 
I = 2,m = 2 symmetric mode. Note that modes with m > 3 
have I > 3, giving them frequencies too large to satisfy the 
instability requirements. The only unstable modes are m = 1 
or m = 2 (see also Watson 1981 for a different argument as to 
why only low m modes are unstable). 

3.3. Shear instability in the shallow water model 

In this section, we describe the results of our stability calcu- 
lation for F > 0. Before presenting the unstable eigenmodes, 
we discuss the general nature of the instability by deriving 
a generalization of Rayleigh's criterion for the shallow water 
model, and the spectrum of r-modes. 

3.3.1. Rayleigh 's criterion 

To derive the equivalent of Rayleigh's criterion for the shal- 
low water equations, we multiply equation Jl 8i by ?/>* /(Q — 
Af2), and integrate over /i from ^ = -1 to fi = +1. Taking the 



imaginary part of the resulting equation gives 
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A non-zero growth rate a requires that the quantity in the 
square brackets change sign somewhere in the domain. Us- 
ing equation ©, we can write this factor as 



, [iFQ r An d ( Q r 

sLh = H — — 

2H dn\H 
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which is the potential vorticity gradient. We see that as sug- 
gested by Dikpati & Gilman (1999), in shallow water it is the 
potential vorticity gradient that must change sign in the do- 
main, rather than the vorticity gradient, thereby including the 
extra physics of vortex stretching in Rayleigh's criterion. 

We evaluate the critical amount of differential rotation s cr i t 
needed for instability by setting (d /dfj,)(fl r /H) = at ^ = 1, 
because by inspection we have found that the critical level 
always appears at the pole as s is increased beyond its critical 
value for our assumed rotation law. Using equations and 
(|6jl for Aft and H, we find that s c „x is given by 



8-40s crit 

■ S crit(3 + S c rit) 



(24) 



or s cth w 8/3F for large F > 40/3. 

We plot i CI it as a function of F in Figure [2] For F w 100, 
only 2% differential rotation between pole and equator is 
required for instability to be possible, an order of magni- 
tude smaller than the 2D case. The dotted curve in Fig- 
ure [2] shows the critical s needed for instability for the ro- 
tation law 1 - ( 1 ?/2)(/i 2 + /i 4 ). For this rotation law, we find 
F = 48(1 -7s cr jt)/(s cr it(19 + 1 lscrit)X giving a ver Y similar crit- 
ical value in the F>1 limit, s cr \t ~ 48/19F = 2.53 /F as com- 
pared to 8/3F = 2. 67 jF previously. In fact, a general rotation 
law An = -s^ n has s crit = (2/F)(n+2)/(n+ 1) for F > 1. The 
critical value s a i t F translates into a critical height contrast be- 
tween equator and pole, using equation of AH /H > 1/3. 
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FIG. 5. — (a) Height contours for the unstable mode shown in Figure^ which has ui = -0.0267 (f e ff = 3.2). The height perturbation is strongly sheared at the 
critical level /x cl j t = (oj/s) 1 / 2 , approximately 20 degrees from the pole. Solid contours are positive, dotted contours are negative, (b) For comparison, the m= 1, 
F = 100, v — 5 r-mode for a rigidly rotating background. The mode frequency is u> = —0.021 1. 





FIG. 6. — As Figure|5] but projected onto a sphere. Colored contours indicate the height perturbation. 



For the more general rotation law AO oc ji n , this becomes 
AH/H > l/(n + l). 

3.3.2. The r-mode spectrum 

We next discuss the r-mode spectrum for rigid rotation 
(s = 0). The effect of increasing F is to reduce the rotat- 
ing frame r-mode frequency (see also the calculations of Lee 
2004). Figure |3 shows the r-mode frequency in the rotating 
frame, 1 + u>, as a function of F for different m and / val- 
ues. On the left, the modes have I = 2,3,4 and 5 from top 
to bottom. On the right, for large F, modes with the same 
l-m have almost the same frequencies. In fact, the r-modes 
of a thin, spherical shell were analysed in detail by Longuet- 
Higgins (1968), who derived an asymptotic form for the mode 
frequency in the limit F ^> 1, 



2mVL 



(2v+\)F l l 2 



(25) 



(eq. [8.33] of Longuet-Higgins 1968), where v = l—m\& anew 
quantum number for large F. This compares well with our 
calculated mode frequencies (also, compare our Figure[3]with 
Figures 2 and 3 of Longuet-Higgins 1968). Later, we will use 
the quantity v e ff = {(2/F 1 I 2 uj)- l)/2 as a label for the modes. 
Using the definition of F, the frequency in the rotating frame 
in the limit F>1 can be written as mui = —m(y/gH /R)/(2v+ 
1), a direct multiple of the surface gravity wave frequency in 
a non-rotating star. 



The eigenfunctions for F > are known as Hough func- 
tions (Longuet-Higgins 1968), shown in Figure |4] for anti- 
symmetric modes with v = 1,3,5 and 7 and F = 100. Note 
that the number of latitudinal nodes in the height perturbation 
h is v - 1 . Modes with v odd are antisymmetric in tp about 
the equator, or symmetric in h. Modes with v even are sym- 
metric in if> about the equator, or antisymmetric in h. Note 
however, that in the geostrophic approximation, h vanishes at 
the equator even for odd v modes (such as those shown in 
Fig-EJ because the Coriolis force vanishes there. This is not 
the case for solutions to the full shallow water equations (see 
Appendix). 

Analytic solutions to the full shallow water equations in the 
limit F» 1 were derived by Longuet-Higgins (1968). In 
the Appendix, we derive a similar solution under the quasi- 
geostrophic assumption. In that case, the perturbation equa- 
tion reduces to the equation for a simple harmonic oscillator 
in quantum mechanics. The eigenfunction is able to propa- 
gate from the equator to a turning point at /i = \/2v+ 1 /F l l A = 
(2v+\)q~ l l 2 , and then evanesces towards the pole. The maxi- 
mum height displacement is close to the turning point. There- 
fore, while rapid rotation (large F) concentrates the eigen- 
function towards the equator, this is less so for large v values. 
For example, the v = 3 mode with F = 100 has a turning point 
at 33° from the pole, and maximum amplitude at 45° (Fig.|4j. 
Modes with v > {F x l 2 — 1)/2 can propagate all the way to the 
poles (v>5forF= 100; v > 8.2 forF = 300). The r-modes 
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FIG. 7. — An example of an unstable mode, the m = 1 antisymmetric mode 
for s — 0.03, F = 100. The solid and dotted curves show the real and imaginary 
part of the height perturbation. The mode frequency is ui = -0.0267, and 
growth time 1.87 s. The critical level for this mode is shown by the vertical 
dashed line (/x CI i t = (lo/s) 1 !' 1 , approximately 20 degrees from the pole). The 
short dashed line shows the potential vorticity profile (normalized so that it's 
maximum value is unity). For this mode, Ugg = ({2/uiF x l-)- l)/2 = 3.2. 



are much less confined than the g-modes, which have a turn- 
ing point at fi = 1/q (e.g. Bildsten et al. 1996). 

3.3.3. Unstable modes 

The discussion so far suggests that shear instabilities will be 
important for the ~ 1 % levels of differential rotation expected 
during Type I X-ray bursts. First, we have shown that a much 
smaller level of differential rotation satisfies Rayleigh's crite- 
rion for large F than for the 2D case. Secondly, the r-mode 
frequencies are much closer to the spin frequency than in the 
2D case, so that it is possible for them to develop a critical 
level at such low levels of differential rotation. The dotted 
curve in Figure [3] shows the critical value of differential ro- 
tation 5 cl ;t needed for Rayleigh's criterion to be satisfied. For 
F <C 1, only the I = 2 mode lies above the dotted curve. We 
saw in SI3.2l that this mode (either with m = 1 or m = 2) is as- 
sociated with the development of dynamical instability as it 
crosses into the corotation region. For the values of F > 100 
expected during Type I X-ray bursts, several modes have fre- 
quencies above the dotted curve, suggesting that dynamical 
instability is possible for large enough differential rotation 
that one of these modes develops a critical level 2 . 

In fact, we do find unstable modes in this region. Figures 
[5] and [7\ show an example of an unstable m = 1 eigenmode 
for F = 100 and s = 0.03. This mode has a rotating frame 
frequency of -0.0267 (which gives f e ff = 3.2), and e-folding 
time 1.87 s. In Figure^ we plot the real and imaginary parts 
of the height perturbation. The vertical dashed line shows the 
critical level, given by — u) = sfi 2 , in this case at 20° from the 
pole. The short-dashed curve shows the background potential 
vorticity profile. The mode frequency is such that the max- 
imum in the potential vorticity profile (ie. the latitude where 
d(fl r / H)/d/j, = Q) occurs close to (but not at) the critical level. 

2 In Figure|5]we are comparing the r-mode frequency based on rigid rota- 
tion (s = 0) with the critical value of s. This is reasonable because the r-mode 
frequencies do not change substantially as s is increased from zero to the 
critical value (e.g. see Fig.lSl. 



FIG. 8. — Mode frequency and growth rate for a variety of modes with 
F = 100. The modes are: m = 1 antisymmetric (solid curves), m = 1 sym- 
metric (dotted), m = 2 symmetric (short dashed), m = 2 antisymmetric (long 
dashed), and m = 3 symmetric (dot-dashed). The mode has a corotation point 
to the right of the short dashed straight line. Rayleigh's criterion is satisfied 
to t he ri ght of the vertical long-dashed line (i cl j t = 0.0234 for F = 100, from 
eq. 1241 ). Note that the mode frequency does not include a factor of m. The 
growth rate does include the factor of m, and is given in units appropriate for 
a star with = 300 Hz (it may be rescaled using a cx f s ). 



This is a common feature of the unstable modes. The max- 
imum amplitude of the unstable mode occurs polewards of 
the critical level, as found for 2D instabilities (Charbonneau 
et al. 1999). Figure |5] shows the perturbation height contours 
for the unstable mode and an r-mode for comparison. The 
height contours are strongly sheared at the critical level. As 
discussed by Dikpati & Gilman (1999), this acts to transport 
angular momentum from the equator to the poles. 

Figure [3] shows the mode frequencies and growth rates for 
F = 100 as a function of s. From equation ([7), we see that 
this value of F is appropriate for the early stages of a burst 
on a star with f s = 300 Hz. To the right of the vertical long- 
dashed line (s > Scrit), the background potential vorticity pro- 
file satisfies Rayleigh's criterion. To the right of the sloping 
short-dashed line, the modes have a critical level. At a given 
value of differential rotation larger than s cl ; t , we find several 
unstable modes. In Figure [8] we start from either s = or 
s = 0.05, and work across the diagram, approaching the coro- 
tation boundary from either the left or right. At some point 
close to the corotation boundary, our code jumps across the 
boundary, moving from either an unstable mode to a stable 
one, or vice versa. Unlike the 2D case shown in Figure^ we 
have not been able to track a given mode across the boundary 
of the corotation region. We expect this is because of numer- 
ical difficulties as the growth rate approaches zero, and the 
critical level becomes singular. A more detailed study along 
the lines of Watts et al. (2003) for the 2D case is needed to 
understand the behavior of the modes close to and at the coro- 
tation boundary. Nevertheless, we can see clearly in Figure 
[8]that there is a correspondence between r-modes outside the 
corotation region and unstable modes inside. 

Figure[8]shows that as s increases, the most unstable mode 
has larger m, larger growth rate, and the mode frequency in the 
rotating frame increases. For s near the threshold, s rj 0.03, 
the most unstable mode is antisymmetric and m = 1. However, 
for s « 0.04, the most unstable mode is symmetric and m = 2. 
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FIG. 9. — As Figure|8] but for F = 300. In this plot for clarity we do 
not distinguish between symmetric and antisymmetric modes, and we do not 
show the stable r-modes (which are located to the left of the short-dashed 
straight line). Alternating peaks in the growth rate as s increases correspond 
to alternating symmetry. We show m = 1 (solid), m = 2 (dotted), m = 3 (short- 
dashed), and m = 4 (dot-dashed) modes. For F = 300, Rayleigh's criterion is 
satisfied for s > Serit = 0.00849, indicated by the vertical long-dashed line. 



FIG. 10. — Similar to Figure|8] but now for a fixed value of s = 0.03, the 
variation in mode frequency and growth rate with F. We show the m = 1 
antisymmetric (solid), m = 2 antisymmetric (dotted), m = 3 symmetric (short- 
dashed), and m = 4 symmetric (long-dashed) modes. For clarity, we show 
either antisymmetric or symmetric for each value of m, choosing the sym- 
metry which gives the largest growth rate in this range of F. The horizontal 
short-dashed line shows the mode frequency below which there is a critical 
level. The long dashed line shows the critical value of F needed to satisfy 
Rayleigh's criterion at * = 0.03 (F cr j t = 75). 



The growth rates are of order ~ 0.00117, or ~ 100 rotation 
periods. For a star with f s = 300 Hz, this is a ~ 1 s _I . Figure|9] 
shows the unstable modes for F = 300. As would be expected 
from Figure [3] more unstable modes are found at this larger 
value of F, and we find instability for a smaller value of s. 
The most unstable modes in this case are those with m = 3 or 
m = 4, and they have larger growth rates than the F = 100 case. 
The curves in Figure [K)] terminate suddenly in many cases 
because we are not able to follow that particular mode any 
further, usually because the code switches to another unstable 
mode. 

Figure ITOl shows the variation with F at a fixed value of 
s = 0.03. The vertical dashed line shows the critical value of 
F needed to satisfy Rayleigh's criterion; below the horizontal 
dashed line a mode has a critical level. As F increases, the 
most unstable mode has larger m. We can understand this by 
inspecting the eigenfunction as m increases. With increasing 
F or s, the maximum in the potential vorticity moves towards 
the equator. Increasing m provides a way to move the maxi- 
mum of the eigenfunction in the same direction. To see this, 
note that the term involving m in equation dl 81 (hidden inside 
the operator C) is most important at the poles, where 1 -fi 2 di- 
verges. In a WKB sense, the effect is to increase the k<j> com- 
ponent of the horizontal wave vector C ~ k\ = Ic^ + ki. The 
mode amplitude increases more slowly away from the pole, 
and the maximum in the mode amplitude occurs closer to the 
equator. 

Shear instabilities in the shallow water model were cal- 
culated by Dikpati & Gilman (1999) for application to the 
Sun, using the full shallow water equations. Our results 
agree well with those of Dikpati & Gilman (1999) for the 
values of s > 0.15 they considered, and help to understand 
some of their results. In particular, they noted that the per- 
turbed fluid velocities were in geostrophic balance except near 
the equator; this provides an additional justification for our 
ageostrophic approximation. They found that the m = 1 and 
m = 2 growth rates decreased for F > 400 (or for their param- 



eter G = 4/F < 0.01), which is roughly where we find that the 
most unstable mode has m > 2. Higher m modes perhaps play 
a role in the solar case. 

4. DISCUSSION 

We have investigated the hydrodynamic instability of lat- 
itudinal shear during the cooling phase of Type I X-ray 
bursts, using a shallow water model. We have shown that 
the geostrophic approximation allows a simple understanding 
of the r-mode spectrum and criterion for instability, although 
at the expense of some accuracy in the eigenfunctions near 
the rotational equator. Considering the simplest shear profile 
symmetric about the equator, we find that a latitudinal shear of 
> 2%, with the equator rotating more rapidly than the poles, is 
linearly unstable to hydrodynamic shear instabilities. Instabil- 
ity is related to low frequency buoyancy-driven r-modes that 
develop a critical layer at which the mode frequency matches 
the local rotation rate in the differentially rotating shell. At 
the same time, a modified version of Rayleigh's criterion, 
based on the potential vorticity, must be satisfied. The un- 
stable modes have low azimuthal wavenumber m, e-folding 
times of roughly a second, and phase speeds close to, but less 
than, the equatorial rotation rate. 

Could the unstable modes be the source of burst oscilla- 
tions? Identification of the burst oscillations with this insta- 
bility would be exciting given the importance of understand- 
ing angular momentum transport in stratified regions in stellar 
interiors. The burst oscillations would then give a remark- 
able opportunity to study angular momentum transport. The 
requirements for the differential rotation profile are that the 
equator should be rotating more rapidly than the poles, by 
about 2%. This is very close to the zonal shear during the 
burst, estimated in |2]to be w 1%. The requirement for insta- 
bility Scrit > 8/3.F for F ^> 1 implies a temperature contrast 
across the star giving AH/H > 0.3. Whether this level of 
temperature contrast is achieved during the burst rise can be 
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addressed by studies of the ignition and propagation of the ini- 
tial burning front. The fact that the estimated level of differen- 
tial rotation is close to the stability boundary may explain why 
burst oscillations are only observed in some bursts. We have 
considered only one particular differential rotation profile, the 
simplest possible profile symmetric about the equator. Differ- 
ent profiles may change the stability boundaries and proper- 
ties of the unstable modes. For example, if ignition occurs off 
the equator, the differential rotation profile may not be sym- 
metric about the equator. It would be interesting to explore the 
effect of the ignition site and spreading on the shear profile. 

The most important question to answer is why only the 
m = 1 mode is excited. Our results suggest that this requires 
some fine-tuning of the level of differential rotation, and is 
only possible for spin frequencies close to 300 Hz. The rea- 
son is that we have shown that unlike the two-dimensional 
case, in which only 1 = 2, m=\ or m = 2 modes are unsta- 
ble (Watson 1981), the shallow water model generally has a 
range of modes that are unstable. The crucial parameter is 
F = (2ilR) 2 1 gH (eq. Q). At the beginning of a burst on a 
star with rotation frequency f s = 300 Hz, F w 100. In this 
case, and for a differential rotation of s 3%, just over the 
threshold for instability, we found that the m = 1 mode is most 
unstable. This changes to the m = 2 mode on increasing the 
differential rotation to s 4%. For larger values of differential 
rotation, or for larger F, the most unstable modes have m > 2. 
This is a problem for the sources with f s w 600 Hz, for which 
F is 4 times larger (F oc / s 2 ). It may be that the non-linear 
evolution of the instability preferentially drives the m = 1 am- 
plitude, and in addition Heyl (2004) points out that the higher 
m modes are much less visible when averaged over the stel- 
lar surface. Another possibility is that the 600 Hz sources 
have m = 2 modes preferentially excited, and therefore have 
f s Rj 300 Hz. However, this may require some fine-tuning of 
the shear profile, especially so that the excited mode always 
has the same value of m in a given source. 

The change in the nature of the instability with spin fre- 
quency is interesting, because it could have implications for 
the spin distribution of the neutron stars in low mass X-ray 
binaries. Chakrabarty et al. (2003) argue that there is an up- 
per limit to the f s distribution of « 700 Hz, consistent with 
observed millisecond radio pulsar spin frequencies, and indi- 
cating that some mechanism, perhaps gravitational radiation, 
prevents accretion from spinning up the neutron star further. 
If the burst oscillations are connected with shear instabilities, 
the lack of high frequency oscillations could be because the 
unstable modes have high m for large f s , and so are unobserv- 
able. However, at this point we do not have a clean prediction 
for the observable spin range. 

The lowest burst oscillation frequency is also interesting 
to consider. EXO 0748-676 has a detected burst oscilla- 
tion at 45 Hz (Villarreal & Strohmayer 2004), significantly 
lower than other sources (the next highest burst oscillation 
frequency is 270 Hz). For a rotation rate this low, the Rossby 
deformation radius is comparable to the radius (F ~ 1), and 
so a large amplitude signal is not expected during the burst 
rise because the rotationally-confined hotspot covers most of 
the star (Spitkovsky et al. 2002). If the observed oscillation is 
due to an r-mode, the actual spin frequency of the star could 
be as much as w 20% larger than the observed frequency in 
EXO 0748-676, since the unstable mode at low F has 1=1. 

We found that the modes driven unstable by the latitudinal 
shear have v > 1 . This makes them more consistent with ob- 
served frequency drifts. Heyl (2004) pointed out that although 



buoyancy-driven r-modes have the required properties to ex- 
plain burst oscillations (low m, a small number of latitudinal 
nodes, they travel backwards relative to the spin), a problem is 
that the surface wave (ie. lowest radial order) associated with 
the burning layer has a frequency that is too large to explain 
the observed drifts. We can see the difficulty from equation 
d25i for the r-mode frequency. In the rotating frame, the fre- 
quency of the r-mode can be written 



mtu = — 




(26) 

where we scale to v = 1, the mode with the fewest latitudinal 
nodes. A factor of 3-5 change in temperature as the burst 
cools gives a frequency drift of Af sa 10 Hz, much larger than 
observed drifts. A reduced vertical lengthscale H » 100 cm 
rather than w 1000 cm is one possible way to reduce the v = 1 
frequency to < 10 Hz, but this requires a mode with more 
than one radial node. Rescaling to m = 1 and v = 3, which 
represents the unstable mode for F = 100 and s = 0.03, and 
using equation (|8) for H, we find 

/ 7m \ 1 /T 9 \ 1/2 
into = — 6.5 Hz — — (27) 

If we assume that the material burns to /x m « 2 rapidly at the 
start of the burst, then the frequency drift expected if the layer 
cools from an initial temperature of Tg = 2 (as limited by ra- 
diation pressure, see Cumming & Bildsten 2000, eq. [12]) to 
Tg = 0.5 is from 6.5 Hz to 3.25 Hz, in better agreement with 
observed drifts. This point is complicated by the fact that the 
frequency drift is likely related to the non-linear development 
of the instability, since it acts to transport angular momentum 
within the burning shell, reducing the background shear. The 
non-linear evolution should be straightforward to follow, by 
direct numerical integration of the shallow water equations. 

Piro & Bildsten (2005) have recently calculated the spec- 
trum of radial modes for detailed models of the cooling at- 
mosphere following a Type I X-ray burst. They find that the 
surface wave associated with the burning layer undergoes an 
avoided crossing with the crustal interface wave as the burst 
cools. The crustal mode has a temperature-independent fre- 
quency, and so this transition terminates the frequency drift, 
bringing the overall drift into agreement with observed values 
for v = 1 . In addition, this gives a well-defined asymptotic 
frequency, resolving the question of why the asymptotic fre- 
quency observed in bursts should be so stable if it is a mode 
frequency (say ss 3 Hz below the spin frequency) rather than 
the spin frequency itself. Calculations of latitudinal shear 
instabilities including the full radial structure of the burning 
layer should be carried out. 

With increasing v, the r-modes should become less visible 
when the flux is integrated over the stellar surface. Many au- 
thors have calculated the lightcurves of rotating neutron stars 
with hotspots of different sizes and locations, including ef- 
fects of gravitational light bending (for direct applications to 
burst oscillations, see Weinberg, Miller, & Lamb 2001; Nath, 
Strohmayer, & Swank 2002; Muno, Ozel, & Chakrabarty 
2002; Bhattacharyya et al. 2005). Heyl (2004, 2005) and Lee 
& Strohmayer (2005) calculate the expected lightcurves for 
low order buoyancy-driven r-modes, but for lowest latitudi- 
nal order modes only. A calculation of the observable ampli- 
tude for the shear unstable modes is needed, especially since 
these modes have a quite different latitudinal variation to the 
r-modes, with maximum height perturbation close to the pole 
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rather than the equator. Interestingly, Muno et al. (2002b) 
find that a hot spot located close to the pole helps to reduce 
the harmonic content of the signal. 

Oscillations were observed during the superburst in 
4U 1636-54 (Strohmayer & Markwardt 2002). Interestingly, 
the value of F will be similar for this case. The superburst ig- 
nition is believed to occur at depths of w 10 12 g cm" 2 (Cum- 
ming & Bildsten 2001; Strohmayer & Brown 2002), where 
the scale height is 50 m. However, taking typical values 
of Ttj = 3 and Ep — 3 MeV, the effective lengthscale for the 
surface wave at these depths is w HiksT /Ep) « 5 m, compa- 
rable to the scale height in a normal X-ray burst. Therefore, 
shear instabilities may be relevant for the case of superbursts 
also. However, the time for burning to spread around the star 
at such high densities is likely to be much shorter than usual 
Type I X-ray bursts, so the level of latitudinal differential ro- 
tation may be small. In this case, another excitation mecha- 
nism for the r-mode needs to be found (such as energy input 
from nuclear burning, also possibly relevant for normal Type 
I bursts, see McDermott & Taam 1987; Strohmayer & Lee 
1996; Lee 2004; Piro & Bildsten 2004). 

Finally, other processes may act to distribute angular mo- 
mentum during the cooling phase of the burst. For example, 
a turbulent shear viscosity might act to damp the modes or 
the background differential rotation. We have ignored the 
effects of radial shear, and any effect of radial angular mo- 
mentum transport on the latitudinal shear profile. Cumming 
& Bildsten (2000) addressed this question, and were unable 
to find a robust hydrodynamic angular momentum transport 
mechanism that acted on a timescale short compared to the 
burst duration. However, short wavelength baroclinic insta- 
bilities should be unstable during the burst tail, and may act 
to transport angular momentum vertically. Most importantly, 



they pointed out that even a weak magnetic field could lead 
to significant angular momentum transport, as it is wound up 
and acts back on the flow (Spruit 1999). The diffusive magne- 
torotational instability can also play a role in vertical angular 
momentum transport within the burning layer (Menou 2004). 
The strength of the magnetic field in the accreted layers is 
uncertain for most X-ray bursters, which show no direct evi- 
dence of a magnetic field. One possibility is that the magnetic 
field in the accreted layers is weak due to screening currents 
in deeper layers (Cumming, Zweibel, & Bildsten 2001). Two 
of the burst oscillation sources are also persistent millisecond 
X-ray pulsars (SAX J 1808.4-3658, Chakrabarty et al. 2003; 
XTE J18 14-338, Strohmayer et al. 2003). The burst oscilla- 
tions from these two sources have unusual properties, such as 
very rapid frequency drift (Chakrabarty et al. 2003), and sig- 
nificant harmonic component (Strohmayer et al. 2003). The 
magnetic field apparently has a direct effect on burst oscilla- 
tions in these two sources. It is possible that a weak magnetic 
field in the other burst oscillation sources plays a role, and in- 
deed a weak magnetic field is known to significantly change 
the stability properties of a shearing flow (e.g. Balbus 1995). 
The possible effect of a magnetic field on the shearing flows 
during Type I X-ray bursts should be considered next. 
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APPENDIX 

CALCULATION OF R-MODES USING THE FULL SHALLOW WATER EQUATIONS 

In 33.11 we made the approximation that the perturbed fluid velocities were in geostrophic balance, and divergence-free. In 
fact, this assumption must break down near the equator where the horizontal component of the Coriolis force due to horizontal 
motions vanishes. Here, we calculate the mode spectrum for a rigidly rotating shell using the full shallow water equations, and 
compare with the results from the geostrophic approximation. 

We consider a uniformly rotating shell with angular speed ft, and define the parameter 3 q = 2Q,/w. The numerical solutions 
of Longuet-Higgins (1968) and Dikpati & Gilman (1999) start with equations for a stream function ^> and velocity potential $, 
which measure the curl and divergence of the horizontal velocity. We adopt a shooting method, in which case it is simpler to 
work with equations for the quantities 




(Al) 



and the height perturbation h (Longuet-Higgins 1968; Piro & Bildsten 2004). In terms of these variables, the perturbed momentum 
and continuity equations may be written 





\b = h 


F m 2 


\ dfj, 1 - M 2 / 




W 1-M 2 J 



(A2) 



and 

(l-f/) — =-b(l-q 2 f/)-qumh. (A3) 
dfi 

We solve these equations using a shooting method, after first factoring out the (1 -/i 2 )'™'/ 2 dependence of the solutions. We start 
integ rating a short distance e <C 1 away from the pole, and choose the ratio of h and b to ensure that the right hand side of equation 
jA3> vanishes there. We then look for solutions with definite symmetry about the equator (either h = or b = at fi = 0). 

Figure lATTl shows the height perturbation for a mode with v = 3, F = 300. We compare the exact solution calculated here with 
the analytic asymptotic formula of Longuet-Higgins (1968) and the quasi-geostrophic solution (also, compare with Figure 10 of 

3 We follow the notation of Bildsten et al. (1996). Longuet-Higgins (1968) defines the alternative A = lu/2Q = l/q. 
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FIG. All. — Eigenfunctions for modes with m = 1 and v = 3 in a layer with F = 300. The solid curve is the solution to the full shallow water equations, the 
dotted curve is the analytic approximation for large F, and the dashed curve is the result of the quasi-geostrophic calculation. 



Longuet-Higgins 1968). The asymptotic expression for h is 

2v+l 



h oc 



2m 



H u -i(.v) + 



1 



2^+2 



H v+l (rf) 



(A4) 



where H v (x) is the Hermite polymonial of degree v, and ry = F l l A ^i. There is a good match except for at the equator, as expected. 
The rotating frame mode frequencies are -0.01719 (exact), -0.01734 (quasi-geostrophic), and -2/(2z^+ \)F X I 2 = -0.01650 (an- 
alytic). The good agreement in the eigenfunctions except at the equator gives us confidence in our quasi-geostrophic results for 
instability. Note that near the pole, the analytic approximation becomes inaccurate because there is not sufficient room in this 
case for the eigenfunction to exponentially decay away before it reaches the pole. 

We briefly discuss the relation of our calculations, which average over the vertical structure, to those of BUC96 and Piro & 
Bildsten (2004), which include detailed calculations of radial eigenfunctions. Perturbing the momentum and continuity equations 
and and manipulating them into a single equation for the height perturbation h, we find 



£h = -^rh 

q 2 



with 



C- 



d_ 



1- 



1 - q 2 jJL 2 d[i 



mq(l+q 2 [i 2 ) 



(1-/J 2 )(1-?V) ' (1-^V) 2 



(A5) 



(A6) 



BUC96 and Piro & Bildsten (2004) solve the equation Ch = —Xh, where A measures the horizontal wavevector (A = /(/+ 1) for 
slow rotation). In our case, A = F/q 2 = u 2 R 2 /gH, because in the thin-shell limit only the shallow water surface wave is present. 
In the slow rotation limit, we recover ui 2 = 1(1+ Y)gH /R 2 = k:\gH as expected. 

ANALYTIC EIGENFUNCTIONS IN THE QUASI-GEOSTROPHIC APPROXIMATION 
We next derive approximate analytic eigenfunctions for the quasi-geostrophic r-modes. Equation (16) with Ail = is 



dfi 



1-M : 



■+F/i +q 



ip = 0. 



(Bl) 



where we have defined the parameter q = 2Q,/oj. Now, following Longuet-Higgins (1968), we write r\ = F 1 / 4 /^, which gives 



d 

drj 



drj 



■+(A-77 2 )= 0, 



(B2) 



VFil-y 2 ) 

where A = —qj \[F = -\JgH /Ruj (defined so that A > since oj < for the r-modes). In the limit F ^ 1, we take 1 - ^l 2 = 
1 - J] 2 /V~F w 1, and neglect the m 2 term, to find 

0+(A-77 2 )^ = O. (B3) 
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The solution is 



ip = H v {rf)e 



V/2 



where H v (rf) is a Hermite polynomial, and A = 2v + 1. This gives the expected mode frequencies, since ui = -y/gH /RA. 
The height perturbation in the geostrophic approximation is oc flip (eq. 1171 ). giving 

h oc r]F' 1/4 H l/ (r})e^' 2/2 . 

Using the recurrence relation for Hermite polynomials, H„ +1 (x) = 2xH n (x)-2nH n - 1 (x), this can be rewritten 



h oc F 



-l/4 g V/2 



H v -i('n) + —H v +i(ri) 

2v 



(B4) 



(B5) 



(B6) 



which differs from the approximation to the full shallow water equations derived by Longuet-Higgins (1968) (eq. IA4I 1 by a 
factor of 2v rather than 2v+ 2 in t he denominator of the second term. 

Finally, note that equation JB3I > is the same as the equation describing a simple harmonic oscillator in quantum mechanics. We 
can now understand the behavior of the r-mode eigenfunctions for increasing v. The eigenfunction is able to propagate to the 
classical turning point at r\ = vA, or fj, = V2^+T/F 1 / 4 = (2v+ l)q~ 1/2 , and then evanesces towards the pole. 
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